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Abstract 

An analytical expression for the maximal Lyapunov exponent Ai in generalized Fermi-Pasta- 
Ulam oscillator chains is obtained. The derivation is based on the calculation of modulational 
instability growth rates for some unstable periodic orbits. The result is compared with numerical 
simulations and the agreement is good over a wide range of energy densities e. At very high 
energy density the power law scaling of Ai with e can be also obtained by simple dimensional 
arguments, assuming that the system is ruled by a single time scale. Finally, we argue that for 
repulsive and hard core potentials in one dimension Ai ~ y/e at large e. 

> 1 INTRODUCTION 

Many theoretical and numerical studies have been devoted to the characterization of chaos in high- 
dimensional systems. Nevertheless, several fundamental issues are not understood. In particular, 
the relation between Lyapunov instability analysis and phase space properties like diffusion of orbits, 
relaxation to equilibrium states, spatial development of instabilities remains to be clarified (see 
Ref. J]]] for a review). 

In this paper we present an analytical estimate of the largest Lyapunov exponent Ai of the Fermi- 
Pasta-Ulam (FPU) model based on the study of modulational instabilities of linear waves. The 
FPU model has played for dynamical system theory a similar role to the Ising model in statistical 
mechanics. Let's just quote some major studies and discoveries motivated by the FPU numerical ex- 
periment: the introduction of the concept of soliton ||, the prediction of the transition to large-scale 
chaos by the resonance overlap criterion [||, the use of KAM perturbation theory and Nekhoroshev 
stability estimates 0, the numerical detection of the strong stochasticity threshold [|]. 

Some attempts already exist 0, of computing analytically Ai in the FPU model. Here we 
present a completely different approach, which emphasizes the relevant role played by some unstable 
periodic orbits related to the representation in Fourier modes |§ (preliminary results were already 
presented in Ref. [[Hj). 

In the present contribution, our analysis will be limited to the asymptotic state, where energy is 
equally shared among all normal modes. In Sect. 2, the modulational instability of a plane wave on 
the lattice is discussed, while in Sect. 3 an approximate analytical expression for Ai is obtained from 
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the growth rates of the modulational instability, and is compared with numerical results. Finally, in 
Sect. 4 it is shown that, for sufficiently high energy density e, the dynamics of the system, in the 
phase as well as in the tangent space, is ruled by a single time scale r(e) and some extensions of the 
results to repulsive and hard core systems are presented. 



2 MODULATIONAL INSTABILITY ANALYSIS 



Denoting by u n (t) the position of the n-th atom (n G {1, . . . , N}), the equations of motion for the 
generalized FPU chain read 



u, 



Un+l + Mn-1 ~2u n + (3 



,2p+l 



[U, 



U n -1 



(1) 



where p is an integer with p > 1, the parameter (3 is fixed to 0.1 for comparison with previous results 
and periodic boundary conditions have been adopted. For sake of simplicity, we first report the 
analysis for p = 1 (i.e. for the usual /3-FPU model) and then we generalize the results to any p-value. 
Due to the periodic boundary conditions, the normal modes associated to the linear part of Eq. ([!]) 
are plane waves of the form 

(2) 



where 9 n (t) = qn — out and q = 2irk/N. The dispersion relation is oo 2 (q) = 4(1 + a) sin 2 (g/2), where 
a = 12/50Q sin 2 (g/2) takes into account the nonlinearity [[II]. The modulational instability of such 
a plane wave is investigated by studying the linearized equation associated to the envelope of the 
carrier wave @. Therefore, one introduces infinitesimal perturbations in the amplitude and phase 
and looks for solutions 



u n (t) = 0„[1 + K(t)] e i[9 " (t)+ ^ ] + 0„[1 + b n (t)] e -*[M*)+<Mt)] 
= 20o [1 + Kit)} cos[qn - ut + if) n (t)] 



(3) 



where b n and ip n are assumed to be small in comparison with the parameters of the carrier wave. 
Substituting Eq. @ into the equations of motion and keeping the second derivative (at variance with 
what has been done for Klein-Gordon type equation |12|, |13|), we obtain for the real and imaginary 
part of the secular term e l (i n ~ U}t ) 



uj 2 b n + 2uip n + 'b n = (1 + 2a) [cos q (b n+1 + 6 n _i) - 2b n ] 

- a (b n+1 + 6 n _i - 2b n cos q) - (1 + 2a) sin q (ip n+ i - ip n -i) (4) 

= (1 + 2a) [cos q (i/j n+ i + i/> n -l) ~ 2^n] 

- (1 + 2a) sin q (b n+l - b n -x) + a (ip n+1 + ^n-i - 2ip n cos q) (5) 



-u ipn - 2ub n + ljj Tl 



Assuming further b n = b e l ^ n + c.c. and ip n 
equations for the secular term e l (Q n - nt ) 



ip e l ^ n ^ + c.c. we obtain the two following 



b n 2 + lu 2 + 2(1 + 2a) (cos q cos Q — 1) - 2a(cosQ - cosg) - 2iip [ton + (1 + 2a) sing sin Q] = (6) 



n 2 



2(1 + 2a)(cosgcosQ - 1) + 2a(cosQ - cosg) + 2ib [con + (1 + 2a) sing sin Q] = (7) 



Non trivial solution for (60,^0) are found only if the determinant vanishes, i.e. if the following 
equation is fullfilled: 



(fi + a/) 2 -4(1 + 2a) sin 2 



q + Q 



(Q - uo) 2 - 4(1 + 2a) sin 2 
4a 2 (cosQ — cosg) 2 



q-Q 



(8) 
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This equation admits 4 different solutions once q (wavevector of the unperturbed wave) and Q 
(wavevector of the perturbation) are fixed. If one of the solutions is complex we have an instability of 
one of the modes (q ± Q) with a growth rate equal to the imaginary part of the solution. Therefore, 
one can derive the instability threshold for any initial linear wave, i.e. any wavevector and any 
amplitude. 

For example, for q = 0, we obtain Q = ± sin ((J/2), which proves that the solution is obviously 
stable since the zero-mode, corresponding to translation invariance, is completely decoupled from 
the others. 

For q = 7r, one can easily see that Eq. (|8|) admits two real and two imaginary solutions if and only if 



2<2 . 
COS — > 



1 + a 
l + 3a 



(9) 



Being a = 12/30 2 , = 3/?(20 o ) 2 = 3/3a 2 (where a = 20o), this formula is equivalent to the formula 
found by Sandusky and Page Hl4l . Moreover, it can be easily shown that the first unstable mode, 
associated to the perturbation, corresponds to Q = 27r/N and therefore the critical amplitude a c 
above which the 7r-mode looses stability is: 



- 2 (#) 




3/3 







1/2 



(10) 



Since for the 7r-mode the energy is simply given by E = N(2a 2 + 4/3a 4 ), we finally get the critical 
energy 

7 cos 2 



2N 

9/3 



sin 



N 



3cos 2 (f) 



1 
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For large values of N, we obtain, in agreement with previous approximate estimations 
the expression 



15, ie 



E r 



7T 



3/3 N 



+ 



(12) 



Above this energy threshold, the 7r-mode is 



that turns out to be quite accurate for N > 50 — 60 
therefore unstable and gives rise to a "chaotic breather" , that has a finite life-time and finally disap- 
pears leading to energy equipartition (this complex relaxation process to equipartition is described 
in detail in Ref. 



3 GROWTH RATES AND LYAPUNOV EXPONENTS 

Again for the 7T-mode, from Eq. (H), one can compute the growth rate a(ii,Q) = lm(Q(Q)) from 
Eq. (|) 



<t(tt,Q) 



\ W4(l + a)(l + 2a) cos 2 — + a 2 cos 4 — - 1 - a - (1 + 2a) cos 2 - (13) 

\ y 2 2 2 



It is plotted in Fig. [L] for two different amplitudes. When the amplitude (or energy density) increases, 
the region of instability extends to a larger region of wavevectors, and in particular the most un- 
stable mode Qmax increases reaching an asymptotic value Q ma , x = 2arccos — 4j ~ 0.427T It is 
important to notice that, for sufficiently high energy, the rescaled growth rate a/a(n, Qmax) does 
not depend on the energy density. This suggests that in the high energy limit a unique time scale is 
present in the system; we will discuss in more detail this point in the following. 
One can also derive the lowest unstable mode, whose asymptotic limit at high energy is Q m in = 
2 arccos ~ 0.67T ~ 1.91 (as previously derived in 0). Fig. ^ reports the stability and instability 
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Figure 1: Shape of the growth rate a(ii, Q). The diamonds correspond to an amplitude a 
the triangles to a = 0.5. 
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regions in the plane of the parameters e = E/N and p(Q) = cos 2 (Q/2). In particular, the marginal 
stability line eo, associated to the vanishing of a(ii, Q), and the line Em = piQmax) are reported. For 
the marginal stability line 6q, our results are fully consistent with a previous theoretical estimation 
reported in ||. 

The above results can be generalized to any power p. In particular, limiting the analysis to the 
7r-mode, the instability condition (P) takes the form 



2Q ^ 
cos — > 

2 



1 + a 



1 + (2p + l)a ' 



(14) 



where a is now given by a = f3 



(2p+l)! 



20o sin f) 



>!(p+l)! V^ UkJlii 2) 

above which the 7r-mode is unstable and for large N we have 



One can therefore derive the critical amplitude 



7T 2 p\(p+ 1)! 

2pPN 2 (2p + 1) 



oc N' 



(15) 



It means that (for fixed N) a c is an increasing function of the power of the coupling potential with 
asymptotic limit linip^oo a c = 0.5. Therefore, in the hard potential limit the critical energy density 
for the 7r-mode is finite (e c = 0.5) for any finite iV-value. It should be noticed that this is not in 
contradiction with the integrability of a one- dimensional system of hard rods, because in the present 
case we have also a harmonic contribution at small distances. 

Finally, we observe that, for the generalized FPU- model, in the high energy limit the growth rate 
<7 p (tc, Q) ps y/a v p (Tr, Q), where a p (7r, Q) is independent of the energy. In summary, for high enough 
values of the energy, 



(Tp(7T, Q) oc e l 



2 2(p+l) 



(16) 



We will see in the last section that this scaling law can be also derived from simple dimensional 
arguments. 

Let us now report an analytical estimation of the maximal Lyapunov. As the system is Hamiltonian, 
Pesin's theorem allows us to identify the Kolmogorov-Sinai entropy Sks with the sum of all positive 



Lyapunov exponents. As the Lyapunov spectrum was shown |TB| to be approximately linear at high 
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Figure 2: The solid line corresponds to the marginal stability curve So, the dashed one to £m = ^(Pmax) 
and the symbols to the points analytically found by Poggi and Ruffo 0. The symbols S and U denote 
stable and unstable region, respectively, separated by the solid line. 



energy, one can relate Sks to the maximal Lyapunov exponent Ai, namely 



A N 

Sks = Y, X * = X ^ • ( 17 ) 
i=i 1 



Let us define a new quantity : the instability entropy 



N/2 

S IE (q) = Y,^ 2 ™/ N ) > ( 18 ) 



where the sum is over all positive growth rates [B|. Our crucial physical ansatz is that Sks — Sie{^) 



From this assumption the following expression for the maximal Lyapunov exponent is readily derived: 

Xi = ^S IE (n) . (19) 

Employing the analytical values for a(ir, Q) reported in Eq. flT5|), Ai can be finally computed. Fig. |3] 
attests that the analytical expression (|I9"D is very accurate. In the same figure the data obtained with 
a completely different approach, developed by Casetti, Livi and Pettini (CLP) ||, are also shown. 
The two methods give almost identical results, apart at very low energy, where the CLP findings [H 



is in better agreement with our numerical data ||10|| . In particular, the low energy scaling turns out 
to be Ai ~ e 2 , instead our formula gives Ai ~ e 1 - 5 . 

As a matter of fact, we obtain a very good agreement between our theoretical estimation and the 
numerical results also for p = 2 and 3 (see Fig. [5]). Plotting Xi(E/N) in a log- log scale we observe 
(see Fig. ^) that at high energy an asymptotic linear behaviour is found for p — 1,2 and 3. In the low 
energy limit, however, for any p the same scaling behaviour should be found, as confirmed by the data 
reported in Fig. ^. The two asymptotic linear behaviors (at high and low e) are separated by a knee 
at intermediate e. An estimation of this transition value can be obtained assuming that the linear 
and nonlinear contributions to uj(q) should be of the same order. We obtain a ~ 1 i.e., an energy 
density of the order of 1/(3 (this is equivalent to the estimation given by mode overlap criterion ||]). 
This knee corresponds to a stochasticity threshold pfl which defines the crossing from weak to strong 
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Figure 3: Comparison of the analytical estimate with numerical results for the maximal Lyapunov 
exponent. The solid curve corresponds to our estimation (Eq. |Tj]), the dashed curve to the CLP- 
estimate using Riemannian differential geometry (see Ref. [§]) and the triangles to our numerical 
results (for the /3-FPU, i.e. p = 1). The dotted curve (resp. diamonds) corresponds to the analytical 
estimate (resp. numerical results) for p = 2 and the dash-dotted (resp. squares) to the analytical 
estimate (resp. numerical results) for p = 3. 



chaos and has been called strong stochasticity threshold (SST) (similar SST-transitions have been 
recently identified also for other models of nonlinear chains [^J). However, we have also found that 



it corresponds to a crossover from extended to a more localized state in tangent space [1C 



It is remarkable to note that Chirikov [21J found similarly the maximal Lyapunov exponent of the 
standard map at high energy by averaging over the phase space the maximal eigenvalue associated 
to the main hyperbolic point. It corresponds in our case to averaging the growth rate (|H|) for the 
unstable periodic orbit q = n over the equilibrium equipartition state (where all modes have the 
same weight). A similar approach is known as Toda criterion p2"| , and although it cannot be used 
cLS cL si gnature of chaos, it can give an approximate estimation of Ai. 

In fact, one can understand this average in a better way by recalling that the modes {7r/2},{27r/3},{7r} 
correspond to the simplest unstable periodic orbits and are also the only three one-mode solutions 
of the (3-FPU problem |J. The calculation of the instability entropies of this three modes shows 
that they are extremely close one to another, contrary to the value for other modes. Recently, it has 
been found that very good estimates for Ai can be obtained by applying the CLP method directly to 
any of the above three one-modes solutions [ . This confirms that the dynamics of these one- mode 
solutions is quite peculiar, because they seem to be "optimal" trajectories in the phase-space to 
estimate the maximal Lyapunov exponent. 

The normalized instability entropy 2SiE{q)/N is plotted against q and e in Fig. |]. This should 



be compared with Fig. 2a in Ref. [p4}| , where a similar plot is reported for the Lyapunov exponent 
computed from the true Jacobian matrix, but taking the linear Fourier modes as approximate orbits 
(moreover the case of fixed boundary conditions is considered). At fixed e both these quantities 
initially grow with q (low g's are more stable than high g's), but then a "stability" region (meaning 
that the growth rate of the instability is significantly smaller) is present around q = 2tc/3 (numerically 
estimated as 0.7 rr by Yoshimura in Ref. p4| ). The "ridge" structure observed in |24j reduces in our 
case to two peaks. The location of the first one depends on e and reaches q ~ 0.562 ir for large energy 
density (it is around 0.57T in the region explored by Yoshimura). The second peak is at (N — 2)ir/N, 
and thus goes to n in the N — > oo limit. We remark that our evaluation of the instability takes fully 
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Figure 4: Normalized instability entropy 2Sie/N as a function of q and e. 



into account nonlinear effects, contrary to Yoshimura's, and that it is completely analytical (apart 
from the necessary numerical evaluation of the roots of Eq. (§)). Mainly, we can confirm that a 
"stability" region (in the previous sense) is present around q = 2tt/3, which corresponds to one of 
the exact solutions found in Ref. 0. 



4 HIGH ENERGY LIMIT 

Here, our claim is that in the high energy limit the temporal evolution of system (H) is essentially 
ruled by a single time-scale. This time-scale can be easily derived by a simple dimensional argument 
(see also [|). At high energy density the potential energy per particle scales as S ~ u 2p+2 , where u 
is the typical size of u n , and similarly for the kinetic energy per particle -j| ~ u 2 . Assuming a virial 
of kinetic and potential energies (this is shown to hold for the FPU model in Ref. 
we can introduce a typical time scale 



V K e 

r\j rs_/ — 

N N 2 



r ~ — ~ e 2 + 2( P +i) (20) 
u 

Therefore, the scaling reported in Eq. (plj) is easily recovered and it is natural to expect that the 
same scaling should hold also for Aj. In particular, for p — 1, Ai scales as e 1 ^ (as also found in 
Ref. §). 

A different scaling, with power 2/3, was found in Ref. |26| at smaller energy densities (closer to the 
strong stochasticity threshold (SST), the knee of Fig. 3). A random matrix approximation |27] was 
invoked to explain these results, which is legitimate only if the process is uncorrelated and one is 
approaching an integrable limit, where the Lyapunov exponent vanishes. Both these approximations 
cannot hold exactly in the region where the 2/3 exponent was found, although a fast decrease of the 
Lyapunov exponent at the knee could mimick the vanishing exponent situation. 
For the /5-FPU model (p — 1) in the high energy limit, an identical scaling has been recently found 



by Lepri |28] for the decay-rate of the time autocorrelation functions of the Fourier modes and in 
Ref. |17j for the inverse of the equipartition time. We have also verified that the normalized Lyapunov 
spectrum Aj = Aj/Ai converges to an asymptotic shape (independent of e) for e — » oo. These results 
confirm our assumption that the dynamics of the system is essentially ruled by an unique time-scale 
for sufficiently high energy density. 
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Figure 5: Logarithm of the maximal Lyapunov exponents for the TCL lattice model as a function of 
the logarithm of the energy density e. Circles represent our numerical data obtained with iV = 64 
and adopting a 6-th order symplectic integrator |31[| , with a relative energy conservation ranging 
from 10~ 12 to 10~ 6 . The solid line refers to a slope 0.5. 



It is important to stress that in the limit of hard core potentials (p — > oo) we find Ai ~ e 1 ^ 2 . Moreover, 
also for potentials that are repulsive (i.e. diverge as (x — c)~ q for x — > c with q > 1, or even steeper), 
we argue that the maximal Lyapunov exponent scales as e 1//2 in the high energy density limit. This 
claim is also supported by the numerical results in Ref |2(J, where the author studies three different 
types of potentials, with a repulsive part at short distances: a diatomic Toda lattice, a truncated 
Coulomb potential (TCL) V(x) = + a; — 1] and a Lennard- Jones (LJ) potential. For the TCL 

potential and the LJ potential, small deviations from the 1/2 power are found by Yoshimura, but 
these tend to reduce as the energy density increases, as shown by our simulations of the TCL model 
reported in Fig. ^. In Ref. [^(J also the equipartition times are reported, whose scaling with e at 
high energy density is always consistent with an exponent —1/2. 

Finally, for a gas of hard spheres in a three dimensional box the analogue of the energy density is 
the temperature T. It is well known that for this model the characteristic time is represented by 
the Enskog collision time t&, which for constant density is proportional to T -1 / 2 ]2"9|] . This result 
suggests that the scaling law r ~ e^ 1 ^ 2 should be valid also for a gas of hard spheres at constant 
density even in three dimensions, as confirmed by recent data for Ai |BU . 



5 CONCLUSIONS 

We have shown that the mechanism of modulational instability of a plane wave (for the three ad- 
missible exact one-mode solutions) is strictly connected to Lyapunov instability and its study gives 
quantitative formula for the maximal Lyapunov exponent. This is still somewhat mysterious, because 
this mechanism deals with the short-time behavior, while the Lyapunov instability arises from an 
infinite-time average. One possible explanation of the mystery is that the considered unstable orbits 
perform a good enough sampling of the phase space, which is uniformly visited by a generic orbit in 
the equipartition state. 

Scaling laws for the Lyapunov exponent at high energy density for generalized FPU chains have 
been obtained using the modulational instability approach, but it is also observed that they can be 
easily guessed by simple dimensional arguments. This points to the existence of a universal time 
scale at high energy, independent of the considered quantity (whether it is the Lyapunov exponent 
or a correlation function). The same universality holds for hard core potentials. 
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